#########################################################################
#To make this figure, you must first run "make_full_panel.R" first.######
#########################################################################
rm(list = ls())
library(plm)
library(AER)
library(foreign)
library(readstata13)

setwd("~/Desktop/QJPS Replication")
data_p <- read.csv("all_countries_ideal_points_panel.csv")

### within and between variance ###

with <- lm(cscore ~ as.factor(cow), data = data_p)

within_resid <-resid(with)

within_var <- sum(within_resid^2)/(length(within_resid) - length(unique(data_p$cow)))

total_var <- var(data_p$cscore)

between_var <- total_var - within_var

within_var/total_var

#######################################################################

data_p <- pdata.frame(data_p, index = c("cow","year"))

a2 <- plm(diff(cscore) ~ lag(cscore,1), data = data_p,model="within", effect="twoways")

rho <- 1 + coef(a2)

rho_var <- rho_period <- rep(NA, 100)

for(i in 1:length(rho_period)){
  
  rho_period[i]<-rho^i
  
  ## delta method
  
  rho_var[i] <- (i*coef(a2)^(i-1))^2*vcov(a2)
}

#####################
#Figure 4, left panel
#####################
par(mar=c(3.5,4.5,3,.5))

plot(1:100, rho_period, type="n", axes=F, xlab = "", ylab="")
lines(1:100, rho_period, col="black", lty=1, lwd=2)

axis(1);axis(2,seq(0,1,.2))
mtext(expression(rho),2,cex=1.5,line=2)
mtext("Years",1,line=2)

mtext("Average Rate of Decay",3,cex=2,1)
mtext(expression(t[1/2]==9.11),line=-.75,cex=1.5)

############################################

half_life <- log(.5)/log(1+coef(a2)[1])

examp <- examp1 <- rep(NA,100)

examp[1]<-tapply(data_p$cscore,data_p$year,quantile)[[112]][2]
examp1[1]<-tapply(data_p$cscore,data_p$year,quantile)[[112]][4]


for(i in 2:100){
  examp[i]<-examp[i-1]*(1+coef(a2))
  
  examp1[i]<-examp1[i-1]*(1+coef(a2))
}

#####################
#Figure 4, right panel
#####################
plot(1900:1999, examp, ylim=c(-1.2,0), xlim=c(1900,1999), type="n", xlab="", ylab="", axes=F)
lines(1900:1999,examp,col="black",lty=1,lwd=2)
lines(1900:1999,examp1,col="black",lty=2,lwd=2)
axis(1)
axis(2,seq(-1.2,0,by=.2))

mtext("Convergence of Constitutional Systems",3,cex=2,1)
mtext("Interquartile Range in 1900",line=-.4,cex=1.5)
mtext("CScore",2,line=2.5)











